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ABSTRACT 

We obtain estimates of Sgr A* accretion flow and black hole parameters by fitting polarized sub-mm 
observations with spectra computed using three-dimensional (3D) general relativistic (GR) magneto- 
hydrodynamical (MHD) (GRMHD) simulations. Observations are compiled from averages over many 
epochs from reports in 29 papers for estimating the mean fluxes F,^, linear polarization (LP) fractions, 
circular polarization (CP) fractions, and electric vector position angles (EVPAs). GRMHD simula- 
tions are computed with dimensionless spins a* = 0,0.5,0.7,0.9,0.98 over a 20,000M time interval. 
We perform fully self-consistent GR polarized radiative transfer using our new code to explore the 
effects of spin a,, inclination angle 6, position angle (PA), accretion rate M, and electron temperature 
Te (Te is reported for radius 6M). By fitting the mean sub-mm fluxes and LP/CP fractions, we obtain 
estimates for these model parameters and determine the physical effects that could produce polariza- 
tion signatures. Our best bet model has a* = 0.5, 9 = 75°, PA = 115°, M = 4.6 x lO^^Mgyear'^ and 
Te — 3.1 X 10^° K at 6M. The sub-mm CP is mainly produced by Faraday conversion as modified by 
Faraday rotation, and the emission region size at 230 GHz is consistent with the VLBI size of 37^as. 
Across all spins, model parameters are in the ranges 9 = 42° — 75°, M = (1.4 — 7.0) x lO~*M0year~^, 
and Te = (3 — 4) X lO^'^^K. Polarization is found both to help differentiate models and to introduce 
new observational constraints on the effects of the magnetic field that might not be fit by accretion 
models so-far considered. 

Subject headings: accretion, accretion disks - black hole physics ~ Galaxy: center - radiative transfer 
- relativistic processes — polarization 



1. INTRODUCTION 

The mass of the Gala ctic Center black hole (BH) 
is M ^ 4.5 • lO^M rr. (IGhez et all [20081 : IR.eid et al.1 
120081 iGiUessen et all l2009f) and the spin is un- 
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It resides at a distance d « 8.4 kpc. 
Because of its proximity, it has been observed in many 
wavelengths: 7-rays, X-rays, IR, (sub-)mm, and radio. 
X-ray bremsstrahlung emission originates from hot gas 
at lar ge radii where the BH's gravity b ecomes impor- 
tant (iNaravan. Yi fc Ma hadeva nI 119951 : iNaravan et al.l 
119981 iShcherbakov fc B aganofl~ l2010[ ) and Compton- 
scattered emission originates from close to the horizon 
(JMoscibrodzka ct al. 2009). X-rays at large radii are 
spatially resolved and have be en used to constrain dy- 
nami cal models for this region (jShcherbakov fc Baganofj 
|2010[) . The sub-mm emission is cyclo-synchrotron emis- 
sion originating from close to the BH. Cyclo-synchrotron 
emission is polarized, and both linear and circular 
polarizations have been observed from Sgr A* at several 
sub-mm wavelengths. The accretion flow was recently 
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resolv ed at 230 GHz poeleman et al.l 120081: iFish et al.l 
120 lit ). General relativistic (GR) effects were deemed 
necessary to explain the small size with full width at 
half maximum (FWHM) of 37/ias. Radio emission is 
also produced by cyclo-synchrotron at larger distances 
from the BH. Relativistic frame-dragging is important 
near the BH, so sub-mm polarized observations and the 
Compton-scattered X-rays might help to constrain the 
BH spin. The goal of the present paper is to model the 
sub-mm in the range of 88 GHz to 857 GHz in order to 
estimate the accretion flow and black hole parameters. 

Sgr A* is a variable source with a variability am- 
plitude routinely reaching 30% in sub-mm. A popu- 
lar approach is to flt simultan eous observation s (e.g. 
lYuan. Ouataert fc Narav an 200 4: IBroderick et al.l l2009l 
in particular, the set from iFalcke et al.l (|1998[ ). However, 
in such an approach, one would use a single simultaneous 
set of observations. However, simultaneous observations 
of fluxes, linear polarization (LP), and circular polariza- 
tion (CP) fractions at several frequencies are not avail- 
able. So we consider non-simultaneous statistics of all 
observations at all frequencies and find the mean values 
and standard errors of quantities at each frequency. 

Numerous accretion flow models have been 
applied to the Galactic Center: advection- 
dominated accretion flow (ADAF) (jNaravan fc Yil 
119951 ) , advection-d ominated inflow-outfl ow so- 
lution (AD IOS) (IB landford fc BegelmanI fT999t) . 
j et-ADAF (lYuan. Ma rkoff fc Falckc 2001), jet 
(jMaitra. M arkoff fc Falckd 120091) . and viscous and 
magnetohydrodynamical (MHD) numerical simulations. 
The quasi-analytical models are useful because there 
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is little expense in changing parameters. However, 
they have a large number of free parameters and also 
incorporate many assumption s that a. r e no t justifi - 
able from first principles (Hu ang et alJ I2008L l2009al ). 
which leads to systematic uncertainties in all fits. 
Numerical simulations require fewer inputs and are 
useful for more quantitative modeling of the plasma 
near a rotating BH. General relativistic (GR) MHD 
(GRMHD) simulations (especially three-dimensional 
(3D) simulations), which are run over a sufficiently 
long duration, are still computationally expensive 
and involv e state-of-the-art codes that are still being 
develo ped (iMcKinnev fc BlandfordI 120091: [Fragile et aLl 



200l INoble fc Krolild I2009t IMoscibrodzka et all 120091: 
Penna et al.l l2010rr Yet, these expensive 3D simula- 
tions are required to model the turbulent disk flow, 
because 2D axisymmetric simulations cannot sustain 
turbulence as shown by generalizations of Cowling's 
anti-dynamo theorem fflide fc Palmed I1982D . Given 
their expense, such 3D GRMHD simulatio ns are limited 
to a region relatively close to the BH (jDexter et al.l 
120091 : IMoscibrodzka et al.l 120091 ) . whereas some emission 
and some Faraday rotation might happen far from the 
BH. So we analytically extend the modeled region out 
to 20, OOOAf , perform radiative transfer, and find the 
best fit to the data. The extension to large radius 
allows us to de fine the electr on temperature more 
consistently (Shar ma et al.l 120071 ). We find a posteriori 
(see Appendix |A| that the simulated polarized spectra 
are not overly sensitive to the details of the analytic 
extensions of density and temperature, but may depend 
on the extension of the magnetic field. 

The radiation close to the BH has been mod- 
eled in Newtonian (jYuan. Quataert fc NaravanI 

l2004fl and quasi-Newtonian approximations 

(iGoldston. Quataert fc Igumenshchevll2005l : iChan et "aD 
I2009D . It has been in o deled i n GR assumi n g un- 
polarized (iFuerst fc Wul l200l iDexter et all ]200l 



Dolcncc et al.' 2009f) and polariz ed (|Broderick et al.l 
2009; Shcherb akov fc Huanj [20111 ) light. Fitting the 
total fiux spectrum might not be sufficient to estimate 
the spin, and naturally one expects polarization to 
provide extra o bservational con s traint s. Spin values 
from Oh. = (iBroderick et al.l |2009[) to a* = 0.9 
(IMoscibrodzka et all I2009D have been estimated. 
We neglect Comptonization (Moscibr odzka et al.l 
|2009| ) and radiation from non-thermal electrons 
(iMahadevanI [1991 iOzel. Psaltis fc NaravanI I2OOOI : 
lYuan. Quataert fc NaravanI 12004^ Emissivities 

are calculated i n the synchrot r on a p proximation 
Legg & Westfold JT968t ISazonovl [l96l IPacholczvkl 
1970: Melrose 197l[ ) with an exact thermal electron 
distribution. Discrepancies with the exact cyclo- 
synch r otron emissivities ( Leung. Gammie fc Noblg 
I2OIII : iShcherbakov fc Huand 1201 ID are negligible as 
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estimated in § |5l Exact Faraday rotation and conversion 
expressions are used (Shcherbakov 2008). 

We compare simulated spectra to observed ones at 
many frequen cies simultaneous l y, ext endin g an approach 
pionee red by IBroderick et al.l (|2009l ) and IDexter et al.l 
(|2009() . We compute the average observed spectra, find 
the deviations of the means, and then compare them to 
the average simulated spectra. In the search for the best 



fit models, we are guided by the value of x^/dof , which is 
the normalized sum of squares of normalized deviations. 
Yet, we leave the exploration of the statistical meaning of 
X^/dof to future work. We search the space of all param- 
eters: spin a*, inclination 9, ratio of proton to electron 
temperatures Tp/Tg [Tp/T^ is reported for radius 6Af), 
and accretion rate M to find the minimum y^ models. 

We summarize the radio/sub-mm observations of Sgr 
A* in § H Qur 3D GRMHD simulations are described in 
§ ISI together with the physically-motivated extension to 
large radii, and the electron heating prescription. We 
run simulations for dimensionless spins a^ = a/M = 
0,0.5,0.7,0.9,0.98. The GR polarized radiative transfer 
technique is described in § |5l 

The set of observations we consider consists of the 
spectral energy distribution (SED) within the 88 GHz to 
857 GHz frequency range, linear polarization (LP) frac- 
tions at 88 GHz, 230 GHz, and 349 GHz, and circular 
polarization (GP) fractions at 230 GHz and 349 GHz. In 
§|6|we discuss our results: the best fit models to the ob- 
servations, the importance of various physical effects in 
producing the observed CP and LP and electric vector 
position angle (EVPA), and image size estimates. We 
produce the simulated images of total and polarized in- 
tensities. Discussion in § |7| compares the results to pre- 
vious estimates, emphasizes the significance of polariza- 
tion, notes the sources of error, and outlines prospects 
for future work. In Appendix VK\ we describe a number of 
convergence tests of our GR polarized radiative transfer 
code and the radial extension of the dynamical model. 
Throughout the paper we measure distance and time in 
the units of BH mass M by setting the speed of light c 
and gravitational constant G to unity. 

2. OBSERVATIONS 

Sgr A* is known to be a highly variable source, yet qui- 
escent models of Sgr A* emission are popular and use- 
ful. Unlike the drastic variations of X-ray and NIR fluxes 
(|Baganoff et al.ll2001l:lGenzel et al.ll2003l) . su b-mm fluxes 
do no t vary by more than a factor of 2 — 3 (jZhao et al.l 
120031 ). We compile the set of observed polarized fluxes 
at each frequency, then we find the mean spectrum and 
the errors of the mean fluxes. 

Previously, the observed flux spect r a were compiled 
bv lYuan. Quataert fc NaravanI (|2004[) : IBroderick et al.l 
(|2009| ). However, both papers summarize a lim- 
ited set of observations and concentrate on simultane- 
ously observed fluxes. Sub-mm flux data reported in 
lYuan. Quataert fc NaravanI (|20()4D consist of a short set 



of observations bv iFalcke et al. 

SMA observations bv lZhao et al. 

(12009) ad ds to these t he rest of S MA total flux data 

(jMarrone eralll2006ai ibl. I2007ll2008l ) . So 6 out of at least 



and one set of 
IBroderick et al.l 



29 papers on sub-mm observations of Sgr A* were taken 
into account. We compute an averaged spectrum based 
on 29 papers reporting sub-mm observations of Sgr A*. 
The reported obs ervations vary in covered period 
from s everal hours (An et al.' '2005") to several years 
(jZhao et al. 2003 : Krichbaum et a l. 2006i). We know that 
variati ons of a factor of 2 may happen within several 
hours (|Yusef-Zadeh et ahl 120091 ). whereas variations by 
more than a factor of several are never observed in the 
sub-mm. So, fluxes observed more than a day apart 
are weakly correlated. The issue of autocorrelation in 
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TABLE 1 

Summary of Sgr A* radio/sub-mm observations 



u IGHzl 



Telescopes 



F. [Jy] 



0.683 ±0.032 ISerabvn et al.lll997 
ll99St [Bower et al.lll999al : IAn et al 



IFalcke et al.r 



LP [%] 



CP [%] 



EVPA 



8.45 



VLA 



120051 ) 



0.871 ±0.012 ISerabvn et al.l997; Falcke et al? 
1998: Bower et al. 2002; Herrnstein et al. 2004,; 
jAn et al. 2005; Yusef-Zadeh et al. 2009) 



-Q.26±0.Q6' ' 

]Bower et al.l 
|1999al) 



-0.62±0.26'' 

( Bower et al. 

l2ooa) 



VLBA, VLA 



22.50 



VLBA, VLA 



0^79 ±0.016 iSerabyii et al. 1997; Falek e et al, 
igga iBower et al.l Il999bl ; iHerrnstein et al 
2004 lAn et al.r 120051: ILu et al.l I200r ' 



Yusef-Zadeh et alJl2007l . l2009D 



1.13 5 ± 0.026 UF alcke et al.l 11998; ILo et al.l 
199a IBower et al] 1999b; Herrnstein et al. 
200l I An et al.l '2005; Shen et al. 2005; 
Krichbauin et al. 2006; Yusef-Zadeh et al. 
2007; Lu et al. 2008; Yusef-Zadeh et al. 2009) 



0T20 ± 0.01" 



1999bl: 



et al.l 



Yusef-Zadeh et al 



20071) 



GMVA, VLBA, 
VLA 



0.55 ± 0.2 2" 

( Bowel et al.l 
1999b: 



Yusef-Z adeh et al 
2007) 



BIMA, MPIfR, 
VLBA, VLA, 
Nobeyama, 
NMA, CARMA 



1.841 ± 0.0 80 I IFalcke et al.l 11998 : 

I Krichbaum et al.l 119981: "JBower et al.1 Il999bl: 
'Doeleinan et ay [20^]; Miyazaki et al. 

2004; Shen lit all 120051: [Krichbauin et al. 
, 200&: .Macguart et al.l I2006t ILu et al.. .200& : 



1.42 



± 0.5 ° 

• et al.l 



Yusef-Zadeh et al.||200g|l 



|Owerj 

1999b'; 

MacQu art et all 

.200&1 



1999bl: 



et al.l 



Shen et al.l 



OVRO, 
CSO-JCMT, 
Nobeyama, 
NMA. IRAM 



1.91 ± 0.15 iScrabyn et al. 1997; Falcke et alT 
1998: Miyazaki et al. 2004; Mauerhan et al.l 
120051 : lYusef-Zadeh et al.ll2009f ) 



Mac guart et akl 
2006) 



145 



Nobeyama, 
NMA, IRAM, 
JCMT 



230 



IRAM, JCMT, 
BIMA, SMA, 
OVRO 



2.28 ± 0.26 ^Falcke et al.l 119981: lAitken et a .1 " 
2003; IMivazaki et al.l 120041 : lYusef-Zadeh et a .1 

2001) 

2.64 ± 0.14 jSerabvn et a l."l997; Falcke et alj 
1993 : [Aitken et aU |2000l; Bower et al. 2003, 
bgoa Zhao et al. 2003; Krichbaum et al. 
'2009: Marrone et al. 2006a, 2007, 2008; 
Doeleman et al. 2008; Yusef-Zadeh et al. 2009) 



34g 



SMA, 
JCMT 



CSO, 



674 



CSO, SMA 



3.18 ± 0.12 (Aitken et al. 2000; An et al.' 
2005: Marrone et al. 2006bl . 120071 . 12008 : 

Yusef-Zadeh et al. 2009) 

1 Marronc Tet al.l 120063 7 
2009) 



7.40 ± .66 

JBower et ^ 
2003, 120051 : 

Marron e et all 
2007, 200I) 



-1.2 ± 0.3" 



^^^^^e^alj 



200gl . l201]li T 



111.5 

Jb 



± 5 3 

e^^n 



20031. I2OO5I: 



857 



CSO 



3.29 ± 0.35 

Yusef-Zadeh et al. 

2.87±0.24 ( Serabyn et al.lll997l: IMarrone et al.l 

2008; Yusef-Zadeh et al. 2009) 



120081 : 



6.50 ± 0.61 

^Marrone et al. 
[20'06b, 2007) 



-1.5 ± 0.3~^ 

fMu noz et al.l 
1201]])) 



IMarrone et al ] 

'2007, 2008) 



146.9 ± 2.2 

( Marrone et al. 
2006b, 2007) 



''The uncertainty of the mean of these quantities is given by instrumental errors. 

''The mean LP at 3.5 mm is computed based on lower and upper sidebands in IMacguart et al.l II2006I ). 
systematic error reported therein. 
'^The mean EVPA at 88 GHz is uncertain due to ±180° degeneracy; e.g. the reported EVPA = 80° could as 

timescales will be addressed in future work. We consider 
the following averaging technique to sample the distribu- 
tions of fluxes. First, we define groups of close frequen- 
cies, the frequencies in each group being different by no 
more than several percent from the mean. There are 11 
such groups (see Table [T]). We exclude papers report- 
ing single frequencies far from the mean of each group. 
In parti c ular, t he 94 GHz and 95 GHz observations of 
iLi et all (I2OO8I): IFalcke e t al.' (Tggl) and the 112 GHz ob- 
servations of IBower et a l. (2001) are excluded. A mean 
frequency is ascribed to represent each group. Then, we 
take all reported observations of each polarization type 
(total flux, LP and CP fraction, EVPA) for each group 
and draw the largest sample of fluxes/polarization frac- 
tions, taking observations separated by at least 24 hours. 
When several fluxes are reported over a period of several 
hours (Yusef-Zadeh et al. 200^, we draw one data point 
from the very beginning of the observation, unless a flare 
is reported to occur at that time. Some of the published 
observations have large error bars. Often such data are 
produced by observing in sub-mm with large beam size, 
but light from Sgr A* is blended wi th dust and other 
sources. In particular, SMT data () Yusef-Zadeh et aU 



The error is based on 0.5% 
well be interpreted as —100°. 

|2009( ). early CSO measurement s (ISerabvn et al.lll997t ). 
and early JCMT measurements ([Aitken et al.ll2000l ) may 
have such issues, so we exclude these data from the sam- 
ple. The interferometric observations, especially with 
VLBI, help to reduce the error from otherwise unreliabl e 
observations, e.g. with BIMA array (jBower et al.l l2001"). 
However, some inconsistencies still exist for simultaneous 
observ ations at the same freque ncy with different instru- 
ments (|Yusef-Zadeh et al.ll2009f l. 

After the sample of fluxes, polarization fractions, and 
EVPAs are found for each frequency group, we com- 
pute the mean and the standard error. The sum- 
mary of results is presented in Table [TJ CP frac- 
tions of -1.2% at 230 GHz and -1.5% at 34 9 GHz 
are based on SMA observations by lMunoz et al.l (|2011l ) 
with the reported ±0.3% instrumental error. Note that 
standard errors in our total flux samples are smaller 



than the error b ars of prior observations (Falcke et al 
199S ; lYuan. Q~ataert fc Naravan 2004: Broderick et al 



2009( 1. but still larger compared t o contemporary single- 



observation instrumental errors (jMarrone et al.l l2007| ) . 
That is, we do not incorporate instrumental error in our 
estimates of standard error of the mean fluxes or LP and 
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EVPA at 230 GHz and 349 GHz (even though the instru- 
mental error of LP at 88 GHz is large) . We do not in- 
corporate the source size measurements (jDoeleman et al] 
120081) in calculating x^/dof, but we check that the best 
bet model is not inconsistent with those observations. 
Figure [1] shows a compilation of the mean quantities and 
their Gaussian standard errors. The data are represented 
by both error bars and the interpolated shaded area. A 
red dashed curve on the F^, plot represents the analytic 
approximation Fi, = 0.248z^°'^^ exp(-(i^/1100)^), where 
flux is in Jy and frequency is in GHz. 

3. THREE-DIMENSIONAL GRMHD SIMULATIONS 

Our radiative transfer calculations take the results of 
simulations of accretion flows onto BHs as input . Thes e 
simulations are similar to those in iPenna et al.l (|2010D . 
Below, we review the methodology. 

3.1. Governing Equations 

We simulate radiatively inefficient accretion flows (RT 
AFs) onto rotating BHs using a three-dimensional fully 
general relativistic code (see § 13. 3|) . The BH is de- 
scribed by the Kerr metric. We work with Heaviside- 
Lorentz units. Our five simulations correspond to dif- 
ferent choices of the dimensionless BH spin parameter: 
a* = 0,0.5,0.7,0.9, and 0.98. The self-gravity of the 
RIAF is ignored. 

The RIAF is a magnetized fl uid, so we so l ve th e 
GRMHD equations of motion (jGammie et al.l 120031 ). 
Mass conservation gives: 



V^(pu^) = 0, 



(1) 



where p is the fluid frame rest-mass density, u^ is the 
contravariant 4-velocity, and V^^ is the covariant deriva- 
tive. Energy-momentum conservation gives 

V;.r,^ = 0, (2) 

where the stress energy tensor T^f includes both matter 
and electromagnetic terms, 

Ti; = {p + Mgas +Pgas + b^W'u, + (pgas + feV2)'5^ - 6^6., 

(3) 

where Ugas is the internal energy density and pgas — 
(r — l)ugas is the ideal gas pressure with F — 4/3. 
Models with F = 5/3 sho w minor differences cornpared 
to models with F = 4/3 (iMcKinnev fc Gammid (20041 : 
iMignone fc McKinnevI 120071 ) . The contravariant fluid- 
frame magnetic 4-field is given by b^ and is related to 
the lab-frame 3-field B'^u via hi" = B^hi^/u*, where 
/i^ = m^'Mi, + S^ is a pro jection tensor and Sji is the 
Kronecker delta function (jGammie et al.|[2003D . We of- 
ten employ b below, which is the orthonormal mag- 
netic fi eld vector in a co moving locally flat reference 
frame (jPenna et al.l 120101 ). The magnetic energy den- 
sity (ub) and magnetic pressure (pmag) are then given by 
Umag = Pmag = bi'b^/2 = 6^/2 = b2/2. Notc that the 
angular velocity of the gas is 17 = u'^ /u^. 

Magnetic flux conservation is given by the induction 
equation 



~g{B'v^ -B^v% 



(4) 



but we use a shock-capturing Godunov method that fully 
conserves energy. So, all dissipation from shocks and nu- 
merical diffusivity (e.g. in shear flows or current sheets) 
is fully captured, as required to study RIAFs. 

In Penn a et al.l (|20ia) . we studied both RIAFs and ge- 
ometrically thin radiatively efficient disks. For the later 
case, a cooling term was added to the energy-momentum 
equation ([2]) to describe radiative losses and keep the 
disk thin. The current set of models are all RIAFs, so 
no cooling term is used. Entropy generated by viscous 
or resistive dissipation is advected along with the inflow 
or transported out via convection or in a wind. 

3.2. Physical Models 

The initial mass distri bution is a n ise n- 
tropic equilibrium torus (iChakrabartil Il985al lbl: 
iDe Villiers. Hawlev fc KrohS 120031 ) with pressure 
p = Kqp"^/^ for Kq = 0.009. The torus inner edge is 
at Tin = 20Af and the maximum density and pressure 
are at i?niax = 65M. We initialize the sol ution so that 
p = 1 at the pressure maximum. As in IChakrabartil 
(|1985aD . the angular velocity distributio n of the ini- 
tial torus is a power law, where for the IChakrabartil 
([T985a) g-parameter we choose q — 1.65 (At large radii 
Vt ^ {r/M)^^.). The thickness of the torus at the 
pressure maximum is then \h/r\ ^ 0.3, where 



\h/r\ 



JJJ\e~7:/2\p{r,9,q^)dAg^dt 
f JJp{r,e,(l))dAe4,dt 



(5) 



where dAgif, = ^—gdOdcf) is an area element in the 
9 — (f> plane, and the integral over dt is a time aver- 
age over the period when the disk is in a steady state 
(see §3.6p . A tenuous atmosphere fills the space out- 
side the torus. It has the same polytropic equation of 
state as the torus, p = Kqp^ , with F = 4/3, and an 
initial rest-mass density oi p — 10~^(r/M)"^/^, cor- 
responding to a Bondi-like atmosphere. The torus is 
threaded with three loops of weak, poloidal magnetic 
fleld: the initial gas-to-magnetic pressure ratio is (3 — 

Pgas,max/7^inag,max — iUU, WUCrC Pgas.max anCL Pniag.max 

are the maximum values of the gas and magnetic pres- 
sure in the torus. This approach to nor malizing the ini- 
tial fi e ld is used in many other stu d ies (iGammie et al.l 
20031 IMcK innev fc Gammiel 120041: IMcKinnev * '2006 
McKinnev fc Naravan 200713; iKomissarov &: M cKinne 
2007HPennaet al.ll20ia) . 

Recent GRMHD simulations of thick disks indicate 
that the results for the disk (but not the wind-jet, which 
for us is less important) are roughly independent of t he 
initial field geometry (jMcKinnev fc NaravanI l2007alrbl: 
Beckwith et al.l l2008al but see also IMcKinnev et al.l 
20121) . The magnetic vector potential we use is given 
by 



l0,N OC 



■ sm 



log{r/S) 
Afield/ (27rr) 



[l + 0.02(ranc-0.5)] 



where v^ = u^/u'^ , and g = Det{g^^) is the determinant of 
the metric. No explicit resistivity or viscosity is included. 



(6) 

with al l other Ag in itially zero. This is the same A^ as 
used in lPenna et al.i (2010) . We use Q = (ugas/ugas,max — 
0.2)(r/M)3/4, g^^^ gg^ Q = if either r < S or Q < 0. 
Here Mg,max is the maximum value of the internal en- 
ergy density in the torus. We choose S = 22M and 



Sgr A* Modeling by 3D GRMHD and Polarized Transfer 
Compilation of observations 



Fy [Jy] 




ISO 



100 



8 15 22 43 88 145 231 349 690 

V [GHz] 



Electric vector position angle [deg] 




349 
V [GHz] 



LP [%] 








CP [%] 




-0.5 


"^ 




-1.0 


M 


^---../^^^ 


[\ 


-1.5 




\ 




M 



231 349 

V [GHz] 



8 15 22 43 88 145 231 349 

V [GHz] 



Fig. 1. — Mean observed SEDs of specific flux F,,, linear polarization (LP) fraction, electric vector position angle (EVPA), and circular 
polarization (CP) fraction. The error bars show la standard error of the mean. The dashed line on the Fi, plot represents the analytic 
approximation i^y(Jy) = 0.248i''''*^ exp(— (!^/1100)^) for frequency u in GHz (not the simulated SED). As noted in Table [Tl the error is 
instrumental for CP at high frequencies and LP at 88 GHz, whereas it is computed from a sample of observed quantities for flux, EVPA 
at all frequencies, and LP at high frequencies. 



Afieid/(27rr) = 0.28, which gives initial poloidal loops 
that are roughly isotropic such that they have roughly 
1:1 aspect ratio in the poloidal plane. The form of the 
potential in equation ([B]) ensures that each additional 
field loop bundle has opposite polarity. Perturbations 
are introduced to excite the magneto-rotational insta- 
bility (MRI). The second term on the right- hand- side 
(RHS) of equation [S] is a random perturbation: ranc is 
a random real number generator for the domain to 1. 
Random perturbations are introduced in the initial inter- 
nal energy density in the sain e way, with an amplitude 
of 10%. In lPenna et al.l ()2010f ). it was found that similar 
simulations with perturt)ations of 2% and 10% became 
turbulent at about the same time, the magnetic field en- 
ergy at that time was negligibly different, and there was 
no evidence for significant differences in any quantities 
during inflow equilibrium. 

3.3. Numerical Methods 

We perform simulations using a fully 3D ver- 
sion of HARM that uses a conser vative shock- 
captu rin g Goduno v scheme ( Ga mmie et al.l 

200l IShafee et al.l [20081 iM cKinnev 2006b; 
Noble et all 120061: iMignone fc McKinnev. ,2007; 



Tchekhovskov. McKinnev fc Naravan I 
McKinnev fc BlandfordI I2009D . We 



20071 : 



use horizon- 

enetrating Ker r-Schild coordinates for the Kerr metric 

Gammie et all 120031: IMcKinnev fc Gammid I2004D . 

which avoids any issues with the coordinate singularity 

in Boyer-Lindquist coordinates. The code uses uniform 



internal coordinates {t,x^^\x^^\x^^^) 

physical coordinates (i, r, 6*, (f>). The radial grid mapping 



') mapped to the 



IS 



r{x 



(1)^ 



i?o + exp(a;(i)), 



(7) 



which spans from i?i„ = O.Qrn to i?out ~ 200Af , where 
rn is the radius of the outer event horizon. This just 
ensures the grid never extends inside the inner horizon, 
in which case the equations of motion would no longer 
be hyperbolic. The parameter Rq = 0.3M controls the 
resolution near the horizon. For the outer radial bound- 
ary of the box, absorbing (outflow, no inflow allowed) 
boundary conditions are used. 
The 0-grid mapping is 



n(2) 



0{x'-^))^ Y{2x'-^> ^1) + {1-Y){2x 



(2) 



1) 



1 



(^/2), 

(8) 

where x^^' ranges from to 1 (i.e. no cut-out at the 
poles) and Y = 0.65 is chosen to concentrate grid zones 
toward the equator. Reflecting boundary conditions are 
used at the polar axes. The (/)-grid mapping is given by 
0(a;(^)) = 2-Kx^^\ such that x'^^^ varies from to 1/2 for 
a box with A(J3 = it. Pe riodic bound a ry co nditions are 
used in the (/)-direction. iPenna et al.l ()2010( ) considered 
various A(/) for thin disks and found little difference in 
the results. In all of their tests, A0 > 7\h/r\ and we 
remain above this limit as well. In what follows, spatial 
integrals are renormalized to refer to the full 27r range 
in (j), even if our computational box size is limited in the 
(/)-direction. For the purpose of radiative transfer, we 
combine two identical regions of size A0 = -k preserving 
the orientation to obtain the span of full 27r. 

3.4. Resolution and Spatial Convergence 
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The resolution of the simulations is Nr x Ng x N^ — 
256 X 64 X 32. This is the fiducia l resolution of 
iPenna et all (|2010D . iShafee et all (|2008D found this res- 
olution to be sufficient to obtain convergence compared 
to a similar 512 x 128 x 32 model. In the vertical di- 
rection, we have about 7 grid cells per density scale 
height. Turbulence is powered by the MRI, which is 
seeded by the vertical c omponent of the magnetic field 
(jBalbus fc Hawle-vlll998() . The characteristic length scale 
of the MRI is the wavelength of the fastest growing mode: 



A 



MRI 



9^"^ 
'"i^' 



(9) 



where u^ is the vertical component of the Alfven speed. 
We find that the MRI is well-resolved in the midplane of 
d isk both initially an d in the saturated state. 

IPenna et al.l (J201C0 studied convergence in N^, Ng, and 
Nci, and found that models with Nr = 256 or Nr = 512, 
Ng = 64 or Ng = 128, and N^ = 64 or N^ = 32 behaved 
similar for disks with similar resolution across the disk. 
O ur resolut ion of the MRI and prior convergence testing 
bv lPennaet al.l (|2010( ) for similarly-resolved disks justify 
our choice of grid resolution. It is currently not compu- 
tationally feasible to perform a similar spin parameter 
study at much higher resolutions, and future studies will 
continue t o explore whether such simulations are fu lly 
converged (jHawlev et al.ll2011l:lMcKinnev et al.ll2012[ ). 

3.5. Ceiling Constraints 

During the simulation, the rest-mass density and inter- 
nal energy densities can become low beyond the corona, 
but the code remains accurate and stable for a finite value 
of b'^/p, b'^/ugas, and Ug^^/ p for any given resolution. We 
enforce b'^ / p < 10, 6^/Mgas < 100, and Ugas//0 < 10 by in- 
jecting a sufficient amount of mass or internal energy into 
a fixed zero angular momentum observer (ZAMO) frame 

with 4-velocity u^ = {—a, 0,0,0}, where a = l/-\/— 5** 
is the lapse. 

We have checked that the ceilings are rarely activated 
in the regions of interest of the flow. Figure [2] shows 
the constrained ratios, b'^/p, fe^/wgas, and Ugas/p, as a 
function of 9 at six radii (r = 4, 6, 8, 10, 12, and 14M) for 
the a* = model. The data has been time-averaged over 
the steady state period from t = 14, OOOM to 20, OOOM. 
The ceiling constraints are shown as dashed red lines. 
The solution stays well away from the ceilings. Thus, 
the ceilings are sufficiently high. 

3.6. Approach to Steady State 

We run the simulations from t = OM to t = 20, OOOM. 
The accretion rate, the height- and 0— averaged plasma 
/3, and other disk parameters, fluctuate turbulently about 
their mean values. The simulation reaches a quasi- 
steady state, when the mean parameter value are time- 
independent. Figure [3] shows the accretion rate and 
height- and 0— averaged 1//3 at the event horizon as a 
function of time for all five models. We take the period 
from t = 14, OOOM to t = 20, OOOM to define steady 
state. 

As shown in IPenna et al.l ()2010D , for disk models like 
the one considered, the disk outside the innermost sta- 
ble circular orbit (ISCO) behaves like the a-disk model 
with a ~ 0.1 across disk thicknesses of h/r ~ 0.05 — 0.4. 




e 



Fig. 2. — Ratios of b^/p, b'^/ug^s, and ngas/p versus 8. Black 
curves correspond to different radii in the flow; from top to bottom, 
r = 4, 6, 8, 10, 12, and 14Af . The data is time-averaged over the 
steady state period of the flow, from t = 14, OOOM to 20, OOOM. 
Numerical ceilings constrain the solution to lie below the dashed 
red lines, but we see that the solution does not approach these 
limits. 

This allows one to accurately infer the timescale for 
reaching "infiow equilibrium," corresponding to a quasi- 
steady flow across all quantities, at a given radius. For 
h/r ~ 0.3 by t ~ 15,000M-20, OOOM (the simulation 
runs till 20, OOOAf , but the initial 5, OOOAf are transients 
not necessarily associated with achieving inflow equilib- 
rium for a sim ple viscous disk), we use the results in 
Appendix B of IPenna et al.l (|2010[ ) and find that inflow 
equilibrium is achieved within a radius of r ~ 25M-30A/ 
for models with a* ~ 1 and r '^ 35M for models with 
a ^ 0. Even for a doubling of the viscous timescale, 
inflow equilibrium is achieved by r ^ 20M-25M depend- 
ing upon the BH spin. This motivates using an analyti- 
cal extension of the simulation solution for radii beyond 
r ^ 2bM as described later in § 14.11 

3.7. Evolved Disk Structure 

FigureHlshows matter stream lines as vectors and num- 
ber density n^ as grey scale map. The large scale vortices 
existing on a single time shot (panel a) almost disap- 
pear when averaged over the duration 6, OOOAf (panel b) 
from times 14, OOOM to 20, OOOM. The density is high- 
est in the equatorial plane on average, but deviations are 
present on the instantaneous map. The ISCO does not 
have any special significance: density and internal energy 
density increase through ISCO towards the BH horizon. 

Figure [5] shows magnetic field lines as vectors and 
comoving electromagnetic energy density ex &^ as a 
greyscale map. The structure of magnetic field at early 
times remembers the initial multi-loop field geometry 
(jPenna et al.ll2010l ). but switches at late times to a heli- 
cal magnetic field structure resembling a split-monopole 
in meridional projection. Such switching of magnetic 
field structure suggests that the final helix with projected 
split-monopole is a natural outcome of any vertical fiux 
being dragged into the BH (although the amount of mag- 
netic flux threading the hole and disk may be chosen by 
initial conditions as described in .McKinncv et al. 2012.). 
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Fig. 3. — Accretion rate and height- and (/(—averaged a = 
Pmag/pgas = 1//9 versus time at the event horizon for all five mod- 
els: a» = (solid light cyan), a» = 0.5 (solid dark red), a» = 0.7 
(long-dashed green), a» = 0.9 (short-dashed brown), and a» = 0.98 
(dot-dashed orange). 

The magnetic field structure of a single snapshot (panel 
a) looks similar to the structure of the linear average 
between 14,000Af and 20, DOOM (panel b). The polar 
region of the flow has the strongest magnetic field. The 
magnetic field lines on Figure [S] illustrate only the di- 
rection of the field's poloidal component. The toroidal 
magnetic field is stronger above and below the midplane 
of the disk outside of ISCO. The toroidal field strength is 
comparable to the poloidal field strength inside the ISCO 
and near the disk midplane. 



Velocity (field) and density (color) 




i-c/M 



Fig. 4. — Stream lines of velocity (red vectors) and number 
density ne (greyscale map) for spin a. = 0.9 at </> = in the 
meridional plane(rc as cylindrical radius): single time snapshot at 
t = 14, DOOM on the upper (a) panel and time average between 
t = 14, OOOAf and t = 20, OOOAf on the lower (b) panel. The cor- 
responding calibration bars of n^ are shown on the right. Number 
density is normalized by its maximum value, and the vectors show 
the poloidal velocity direction. 



4. DYNAMICAL MODEL BASED ON SIMULATIONS 
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Fig. 5. — Magnetic field lines (red vectors) and comoving electro- 
magnetic energy density oc b^ (greyscale map) for spin a* = 0.9 at 
</> = in the meridional plane(rc as cylindrical radius): single time 
snapshot at i = 14, OOOM on the upper (a) panel and time average 
between t = 14, OOOM and t = 20, OOOM on the lower (b) panel. 
The corresponding calibration bars of comoving b^ are shown on 
the right. Magnetic field energy density is normalized by its max- 
imum value. The magnetic field lines illustrate only the direction 
of the field's poloidal component. 

We now discuss extensions of the numerical simula- 
tions, which we need to perform radiative transfer com- 
putations. We extend the simulations to large radii and 
define the electron temperature. 

4.1. Extension to Large Radius 

The flow is evolved in a quasi-steady state for 6, OOOAf 
from 14, OOOM until 20, OOOM, which corresponds to 8 
orbits at r = 25M. The flow is not sufficiently settled 
at larger radii. However, outside 25M, some Faraday 
rotation and emission might occur. So, we extend the 
dynamical model to larger radii (i.e. r > 25Af) in a 
continuous way and check (see Appendix \^ how vari- 
ations of our large radius prescriptions change the re- 
sults of radiative transfer. The outer radial boundary 
of radiative transfer is situated at r = 20, OOOM. The 
profiles of number density Ue, internal energy density 
Ugas, magnetic field b, and velocity v are extended as 
power-laws until radius r — 20, OOOM. The power-law 
index for number density /? is obtained by matching the 
known value Ue = IS Ocm"'^ at about 1.5" w 3 • lO^M 
(jBaganoff et al.l l2003| ) and the average ne,cut value at 
r = 25M in the equatorial plane for each model. The 
value of (3 may be different for different models. The 
radial flow velocity Vr is then obtained from the continu- 
ity relation in the equatorial plane n^Vrr'^ ~ const. The 
power-law of internal energy density Ugas is obtained in a 
similar way by matching the values Tg = Tp = 1.5 • 10^ K 
and Uf., = 130cm~^ at distance 3 ■ lO^ M (jBaganoff et al.l 
[2003t IShciierbakov fc Baganofa [20T0I ). The meridional 
physical velocity is extended as Vq cx {r/M)^'^/'^ and 
toroidal velocity as vi cx {r/M)~^^'^ to approximately 
match the power law between 15M and 25M, where 
the relationship uj sa u^^^^/g~ is used to connect the 4- 
velocity components with physical velocity components. 
All components of comoving magnetic field are extended 
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as br,be,b^ ex ir/M) ^, which appears valid across a 
diverse set of GRMHD models (|McKinney et al.ll2012[ ). 
This power-law slope corresponds roughly to equipar- 
tition of niagnetic field energy density, since constant 



fraction niagnetic field is h 



(r/M)-^ for 



n (x (r/M)~^. Exploration of various extensions of the 
magnetic field will be the topic of future studies. 

After defining the extension power-laws for quantities 
in the equatorial plane, we extend the quantities radially 
at arbitrary 9 and (j) in a continuous way. For example, 
for density at arbitrary 9 and (j) ^^nd at r > 25M we have 

n,{r, 9, 0) = ne(25M, 9, 0) (^^'^ "^ , (10) 

where ne{25M,9,(j)) is taken from the simulations. We 
similarly extend other quantities. As shown in Ap- 
pendix |X1 small variations in power-law indices of num- 
ber density and temperature have little influence on ra- 
diation intensities and linear/circular polarization fluxes, 
but variations of magnetic field slope can make a substan- 
tial difference. 



4.2. Electron Temperature 

Neither the proton Tp nor the electron Tg temperature 
is given directly by the simulations. However, it is cru- 
cial to know the electron temperature Tg to determine 
the emission. Our solution is to split the total internal 
energy density Mgas, given by the simulations and their 
power-law extension, between the proton energy and the 
electron energy. The energy balance states 



-''gas ^PrQ 



^e,g 



P 



P 



= CnksTr, + CeksTi 



S^e, 



(11) 



where Cp = 3/2 and Ce > 3/2 are the respective heat 
capacities, p is the rest-mass density, and ks is Boltz- 
mann's constant. The difference of temperatures Tp — Tg 
is influenced by three effects: equilibration by Coulomb 
collisions at large radii, the difference in heating rates fp 
and fe of protons and electrons operating at intermedi- 
ate radii, and the difference in heat capacities operating 
close to th e BH. Radiative coo ling is ignored since, ac- 
cording to iSharma et al.l ()2007[ ). the radiative efficiency 
of the flow is negligible for realistic M < 10~^MQyear~^. 
The relevant effects can be summarized by the equation: 



djTp - T,) 
dr 



-Vc{Tp - Te) + 



(12) 



1 /p 1 /e ^ d(Mgas/p) 



Cp fp + fe c'g fp + fej ksdr 



where 



= 8.9 • 10" 



n 



3-10 



10 



-3/2 



10^ 



(13) 



is the non-relativistic temperature equilibration rate by 
collisions (Slikarofsk y et al. 1966) , all quantities being 
measured in CGS units. We consider protons to always 
have non-relativistic heat capacity and collisions to al- 
ways obey the non-relativistic formula. The magnitudes 
of errors introduced by these simplifications are negligi- 
ble. The exact expressions for total electron heat capac- 



ity and differential heat capacity are approximated as 

_We,g/p 3 0.7 +29e 



ksTe 2 0.7 - 

d{Ue,g/p) _ „ 



0.735 



(14) 
(15) 

(16) 




^ kedT, (0.7 + 9eY 

correspondingly, with the error < 1.3%, where 

a ksTe 
^^ ^ 2 

is the dimens ionless electron tem perature. It was re- 
cently shown (jSharma et al.lf2007l ) that the ratio of heat- 
ing rates in the non-relativistic regime in a disk can be 
approximated as 

(17) 

with coefficient C. This formula is adopted in the rel- 
ativistic r egime as well, since no better prescription is 
available. ISharma et al.l (120071 ) found the value C = 0.33 
in simulations, whereas we find C = 0.36 — 0.42 for the 
best fit models (see Tableland §11]). 

The proton and electron temperatures are determined 
at each point in the following way. We first take a sin- 
gle snapshot of a simulation with spin a* and extend the 
flow quantities to r = 20, OOOAf (see § Hj]). Then we 
compute azimuthal averages of radial velocity w^, num- 
ber density Ug, and Wgas/p at the equatorial plane, ex- 
tend them as power laws to rout = 3-lO^M, and solve the 
equations (|llll2p from rout down to the inner grid cell 
point. Temperatures are set t o T^^ = T^ = 1.5 • 10^ K 
at ro ut (jBaganoff et al.l l2003t iShcherbakov fc Baganofj 
l2010f ). On the next step we compare the values of Wgas/p 
to the calculated Tg and Tp and determine the functional 
dependence T^ = Te{u^a.s/p) and Tp = Tp{ugas/p)- At 
each point of the simulation (including off the equator), 
we draw temperatures from this correspondence. That is, 
GRMHD simulation directly provides Wgas and p at the 
equatorial plane, so the function Tg = T'e(ugas/p) gives 
Te at each point in space. Typical profiles of proton and 
electron temperatures are shown on Figure [6l Tempera- 
tures stay equal until r ~ 10'*M due to collisions, despite 
different heating prescriptions. Within r = 3 • lO'^Af the 
timescale of coUisional equilibration becomes relatively 
long and electrons become relativistic, thus Tg deviates 
down from Tp. The electron and proton temperature 
profiles in the region r < 20,000M are used to conduct 
the radiative transfer. For a given accretion rate, we find 
that there exists a unique dependence of the ratio of tem- 
peratures Tp/Te (measured at r = 6M at the equator) 
upon the heating coefficient C, so that we can use Tp/T^ 
and C interchangeably. 

5. GENERAL RELATIVISTIC POLARIZED RADIATIVE 
TRANSFER 

5.1. Description of Radiative Transfer 

Now we convert the dynamical model of the ac- 
cretion flow into a set of obs ervable quantities us- 
ing polarized radiative trans fer (jBroderick et al.l 120091 : 
Shcherbakov fc Huang] l20ll . We closely follow 

Shc herbak ov & Huang (2011") for the transfer technique. 
Similar to .Huang et al., (2009a^ . we define the polarized 
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proton temperature T^ (red) & electron temperature Tg (blue) 




Fig. 6. — Temperatures of protons Tp (upper red line) and 
electrons Te (lower blue line) for the dynamical model with spin 
a» = 0.5 giving the best fit to polarization observations (see Table|2] 
and §[6)1. 

basis in the picture plane, where one vector points North, 
another vector points East, and the wavevector points 
towards the observer. We parallel transport this basis 
in the direction of the BH and do the radiative trans- 
fer along the ray in the opposite direction towards the 
observer. At each point along the ray we go into the 
locally-flat comoving frame, calculate the angles between 
the magnetic field and basis vectors, and compute the 
Faraday conversion, Faraday rotation, emissivities, and 
absorptivities. 

Radiative transfer involves shooting a uniform grid of 
Pjv X Pn geodesies from the picture plane down to the 
black hole. The total polarized fluxes are computed by 
integration of intensities along each ray backwards to 
the picture plane. We found that Pn = 111 is good 
enough to compute the spectrum (jDexter et al.l 120091 
used Pjq = 150). For radiative transfer we employ all 
3D dat a in each numerical si mulat ion snapshot and, fol- 
lowing iMoscibrodzka et al.l (|2009| ) , perform multilinear 
interpolation in three dimensions for the quantities in 
between the grid points. We make no approximations in 
the use of spatial 3D data. We self-consistently take into 
account the evolution of the numerical simulation as the 
light geodesies travel around the BH. Since it is too time- 
consuming to look up simulation data over a long period 
of time, we only evolve the simulation between t — At and 
i -I- Ai to get a spectrum at time t + 20, OOOM. The off- 
set 20, OOOAf appears, since the picture plane is located 
20, OOOAf away from the BH center. The extension to the 
large radius outside 25 M, however, is not evolved with 
time. It is taken to be that of a single snapshot at time 
t. The snapshot at times t — At and t + At are taken to 
represent the numerical simulations at earlier and later 
times, respectively. We find that At = 6QM is large 
enough to achieve accurate simulated spectra. The total 
fluxes are found at regular time intervals within period 
of quasi-steady accretion from 14, OOOAf tiU 20, OOOAf, 
e.g. for t = 14, OOOAf, 14, 300Af, ..., 19, 700Af, 20, OOOAf. 
We compute iVpcriods = 21 spectra over the quasi-steady 
accretion phase and average them to find the mean sim- 
ulated spectra. To compute the polarized fluxes we take 
the integration domain in the picture plane to be a square 
with a side 



a[Af] = 16 + 2 



600 

i^[GHz] 



1.5 



in the units of Tg = Af , where frequency v is in GHz. This 
square is centered at the BH. The size based on Equa- 
tion (TTSl) is larger than the photon orbit visible diameter 
dph ~ 10.4Af and follows the intrinsic size depende nce 
on frequency (jShen et al." 2005; Doele man et al.ll200 8') at 
low frequencies. An important radiative transfer param- 
eter is the distance from the BH, where intensity integra- 
tion starts. The dependence of synchrotron emissivity 
on temperature and magnetic field strength is so strong 
that it overwhelms the sole effect of gravitational red- 
shift close to the BH. We obtain accurate results in the 
sub-mm for computation out from r^i^ — l.Olrn, where 
rn = Af (1 -|- ^yl — a^) is the horizon radius. To quantify 
the needed accuracy of computations, we define a quan- 
tity xlf/dof in Appendix \X\ We conduct multiple tests 
of radiative transfer convergence for best fit models at 
each spin. In Appendix \^ we justify the chosen values 
of radiative transfer parameters Pn, At, A^poriods, ''min, 
etc. 

Our calculation of plasma response is different from 
iShcherbakov fc Huand (|201lD . They offered a way to find 
exact emissivities, absorptivities, Faraday rotation, and 
conversion coefficients for thermal and other isotropic 
particle distributions. Here, for simplicity, we employ 
fitting formulas for Faraday rotation and Faraday con- 
version and synchrotron approximation for emissivities 
for a thermal plasma. We define 



X 



3 vb^^ sin0B 



(19) 



where 6^ is k-b angle, 7 is electron gamma factor, and 
i^B = ebl{2-K7nec) is the cyclo t ron frequency. Then fol- 
lowing lleg^^&Sitfoii ()1968D : lMelrosi()1971[ ). we write 
down emissivities in the f , Q, and V modes as 



V3e2 

£/ = — vbsixyOb 

2 c ji 

V3e2 . 
eQ = — vbuvhOb 



2 e^ 

£V = —F^ — VB COS ( 

V3 c 
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+00 



d-fN{^)XK^,:,{X) 
0(7 X 



(20) 



p + ca 

XKy,{X) + / dzKy,{z) 
Jx 



(18) 



Here K^ (x) is the Bessel function of the 2nd kind of or- 
der z. We employed lEEE/IAU ciefiniti ons of Stokes Q, 
U, and V (jHamaker fc Bregmanlll996l ). and we define 
counter-clockwise rotation of the electric field as seen 
by the observer as corresponding to positive V^ > - 
as also chosen in'S hcherbakov fc Huand (|2011l) . So, the 
si gn of the V emissivity fEq.J^Ol) is opposite to the sign 
in iRvbicki fc LightmanI (119671). A v a riatio n of emissiv- 
ity for mulas (|19l20p exists: ISazo"novl (|1969( ): IPacholczvkl 
(| 19701 ) define X = 2:// (31/3(7 - l)'^sm9B), integrating 
over particle energy instead of 7. This approximation 
appears to give significantly larger errors at low particle 
energies. 

Next, one needs to identify the accurate thermal par- 
ticle distribution A^(7). Various A^(7) correspond to var- 
ious synchrotron approximation s. The ul trarelativistic 
thermal approximation (|Pacholczvkl I197Q : .Huang et aTl 



10 



Shcherbakov, Penna & McKinney 



l2009al) has the simplest distribution N{'-f) = exp(— (7 — 
l)/^e)(7~ l)^/2/^e- However, the exact thermal distri- 
bution of particles 



iV(7) = 7/? 



-exp(-7/6'e 



Oe.Ko 



") 



(21) 



allows for more precise computation of radiation. Syn- 
chrotron cmissivities based on the equations (|19I20I) 
with the exact thermal distribution (pij) agree with 
the exact cyclo-synchrotron cmi s sivities £j, eg, and ey 
(iLeung. Gammie fc Nobld 120111 : iShcherbakov fc Huang] 
120111) to within 2% for typical dynamical models and fre- 
quencies > 100 GHz. Emissivities integrated over the ul- 
trarelativistic thermal distribution typically have ~ 10% 
error. 

Thermal absorptivities are found from emissivities 
(Eq. [201) via Kirchhoff 's law 



a/,Q,v — £i,Q,v/B^, 



(22) 



where B^, = 2fcsT'e^^/c^ is the source function for low 
photon energies {hv <^ ksTe). Faraday rotation pv 
and Faraday conve rsion pq coefficients are taken from 
IShcherbakovl(l2008l) : 



pv = g{Z) 



PQ^fiZ) 



rrii^cv^ 






(23) 



sin^ e. 



Here 



and 



Z = 



V2 



smf 



103^ 



(24) 



g{Z) = 1 - 0.11 ln(l + 0.035Z), 
/(Z)= 2.011 exp 



4.7 



(25) 



cos 



exp 



Z 



1.2 



2.73 



- 0.011 cxp - 



47.2 



are the fitting formulas for deviations of pv and pq from 
analytic results for finite ratios of vb/^- The deviation of 
f{Z) from 1 is significant for the set of observed frequen- 
cies i^, temperatures 0e, and magnetic fields found in the 
typical models of Sgr A*. These formulas constitute a 
good fit to the exact result for the typical parameters of 
the dynamical model (Shcherbakov 2008). 

Polarized radiative transfer can take much longer to 
perform compared to non-polarized radiative transfer 
when using an explicit integration scheme to evolve the 
Stokes occupation numbers Nq, Njj, and Ny. Large 
Faraday rotation measure and Faraday conversion mea- 
sure lead to oscillations between occupation numbers. 
One of the solutions is to use an implicit integration 
scheme, while another solution is to perform a substi- 
tution of variables. In the simple case of Faraday rota- 
tion leading to interchange of Nq and Njj , our choice of 
variables is the amplitude of oscillations and the phase. 
Thus, the cylindrical polarized coordinates arise as fol- 



lows: 



Nq = Nqu cos ( 
Nu = Nqu sin q 



(26) 



Then, the amplitude Nqu slowly changes along the ray 
and the angle (j) changes linearly, and this translates 
into a speed improvement. In the presence of substan- 
tial Faraday conversion, the polarization vector precesses 
along some axis on a Poincare sphere, adding an inter- 
change of circularly and linearly polarized light. So, po- 
lar polarized coordinates are more suitable in this case: 



iVg^A^poiCOS^sinV', 
Nu — Npoi sin (/) sin ip, 
Nv = Npoi cos ip, 



(27) 



where TVpoi is the total polarized intensity, <f> angle 
changes are mainly due to Faraday rotation, and ip angle 
changes are mainly due to Faraday conversion. The ap- 
plication of this technique speeds up the code enormously 
at low frequencies oi v < 100 GHz. 

5.2. Search for the Best Fits 
We define x^/dof quantities to discriminate between 
models. We define Xf f'^'' fitting total fiuxes as 



xl 



E 



{F^. 



..oh: 



■)' 



a{FY 



(28) 



for the set of 7 frequencies v = 88, 102, 145, 230, 349, 680, 
and 857 GHz, where a{F) are the errors of the means. 
We incorporate LP fractions at 88, 230, and 349 GHz and 
CP fractions at 230 and 349 GHz to obtain 



xl- 



E 

2 

-E 

1=1 



(LPi 



LP,; 



.? 



o-(LP)2 

(GPi^sini ~ GPi^obs) 



(t(CP)2 



(29) 



Then we define dof (as degrees of freedom) to be Aoip = 
7 - 3 = 4 for flux fitting and dof = 12 - 3 = 9 for fit- 
ting all polarized data. The quantity x^/dof would be 
drawn from x^ statistics if ct-s were the true observa- 
tional errors and if the observed fluxes were drawn from 
a Gaussian distribution. However, for the purpose of the 
present work, we only employ x"/dof as a measure of 
fitting the data. That is, lower x^/dof indicates better 
agreement with the data. We do not attempt to ascribe 
any statistical meaning to the quantity x^/dof. 

We explore models with 4 parameters: spin a*, inclina- 
tion angle 6, accretion rate M , and the ratio of proton to 
electron temperature Tp/T^ [Tp/T^ is reported for radius 
r = 6M). For the radiative transfer calculations, the 
density from the simulations is scaled to give the desired 
accretion rate. 

6. RESULTS 

In previous sections, we described our compiled ob- 
servations, GRMHD numerical simulations of the fiow 
structure, our method for obtaining the electron temper- 
ature, and our method for polarized radiative transfer. In 
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TABLE 2 
Properties of the best fit models with different spins. 



Model 


Inclination an- 
gle 0, deg 


Spin position 
angle PA, deg 


Heating con- 
stant C 


Ratio Tp/Te at 
6A/ 


Electron Te at 
6A/, K 


Accretion rate 

M, Mgyr-l 


spin a, = 


42.0 


171.0 


0.42107 


15.98 


3.343 ■ 10"' 


7.005 ■ 10-» 


spin a, =0.5 


74.5 


115.3 


0.37012 


20.14 


3.087 ■ 10"' 


4.594 ■ 10~» 


spin a» =0.7 


64.5 


84.7 


0.37239 


20.16 


3.415 ■ 10"' 


2.694 ■ 10-» 


spin a» = 0.9 


53.5 


123.4 


0.39849 


18.16 


4.055 ■ 10"' 


1.402 ■ 10"** 


spin a* = 0.98 


57.2 


120.3 


0.41343 


17.00 


4.190 ■ 10"' 


1.553 -lO-* 


spin a, = 0.5 
short period 1 


70.0 


79.3 


0.38934 


18.50 


3.334 ■ 10"' 


3.513-10-8 


spin a, = 0.5 
short period 2 


72.8 


113.1 


0.40507 


17.31 


3.541 ■ 10"' 


3.452 ■ 10-s 


spin a, = 0.5 
short period 3 


73.4 


57.4 


0.37302 


19.87 


3.125 ■ 10"' 


3.897- 10-s 


spin a, = 0.5 
short period 4 


74.4 


115.4 


0.36147 


20.95 


2.978 • 10"' 


4.508 ■ 10-s 


spin a, = 0.5 
short period 5 


71.9 


95.7 


0.37420 


19.79 


3.137 ■ 10"' 


5.334 ■ 10-** 


spin a» = 0.5 
short period 6 


76.4 


116.7 


0.38853 


18.59 


3.320 ■ 10"' 


6.080 ■ 10-8 


spin a* = fast 
light 


41.4 


187.5 


0.41929 


16.09 


3.322 ■ 10"' 


7.044 ■ 10-** 


spin a, = 0.5 
fast light 


72.7 


105.9 


0.39804 


17.83 


3.447 ■ 10"' 


3.957-10-** 


spin a, = 0.7 
fast light 


59.4 


131.8 


0.35708 


21.62 


3.204 ■ 10"' 


2.966 ■ 10-** 


spin a, = 0.9 
fast light 


53.3 


123.3 


0.40215 


17.86 


4.116 ■ 10"' 


1.340 ■ 10-** 


spin a* = 0.98 
fast light 


57.7 


119.6 


0.41720 


16.73 


4.246 • 10"' 


1.515-10-** 



Note. — Mean values are shown for ratio Tp/Te, electron temperature Te, and the accretion rate M. These are the simple means over 
all simulation snapshots, which were employed for radiative transfer in a particular model. The values of x^/dof ranges from 2—5 across 
all models. 
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Fig. 7. — Fits to the observed fluxes, LP and CP fractions by best models for each spin. The inclination angle 9, accretion rate Af , ratio 
of temperatures Tp/T^ were adjusted for each spin to minimize x'^/dof . Fits to total flux Fi, are in the upper left panel, LP fraction in the 
lower left, and CP fraction in the lower right. Shown are the best fit models with spin a, = (short-dashed brown), spin a» = 0.5 (solid 
dark red), spin a, = 0.7 (long-dashed green), spin a» = 0.9 (solid light cyan), and spin a» = 0.98 (dot-dashed orange). The upper right 
panel shows the dependence of EVPA on frequency for the best models. Note, that EVPAs are not included into our fitting procedure. The 
thick blue curve represents observations. Simulated EVPA curves are arbitrarily shifted to approximate EVPA at 349 GHz. The addition 
of an external (to the emitting region) Faraday rotation screen helps to fit EVPA(349 GHz) — EVPA(230 GHz). 

this section, we discuss our results for accretion flow and BH parameters, as guided by a minimization of x^/dof 
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for our model applied to the observations. 

Figure [7] shows best fits to observations by models 
with five different spins. Inclination angle 0, accretion 
rate M, and heating coefficient C were adjusted to reach 
the lowest x^/dof. Fits to fluxes F^, (upper left) are 
not substantially different, although models with higher 
spins fit better at high frequencies. Larger deviations 
can be seen on LP (lower left) and CP (lower right) 
plots. Models with high spins require lower accretion 
rate (i.e. density) to fit the fiux spectrum. As a conse- 
quence, they are not subject to Faraday depolarization, 
which leads to a decrease of LP at low v, and the mod- 
els end up having larger linear polarization fractions at 
88 GHz. Not all models reproduce the observed decrease 
of mean LP fraction between 230 GHz and 349 GHz 
groups. The discrepancies in fitting the CP fraction are 
also large: all the lowest x^ models give |GP| < 1.5% 
at 349 GHz. The best bet model with spin a* = re- 
produces LP and CP fractions well, but fails in fitting 
the total flux. Most solutions predict the wrong sign 
of the EVPA(349 GHz) - EVPA(230 GHz) difference, 
which could be fi xed with stronger ma gnetic field (e.g. as 
seen in models bv lMcKinnev et aLll2012) to yield stronger 
Faraday rotation. In sum, crude agreement of simulated 
polarized spectra to the observed ones was achieved, but 
the improved dynamical models may be needed for better 
fits. 

We now isolate the physical effects responsible for the 
observed polarized quantities for our best bet model with 
spin a, = 0.5 that has the lowest x^/doi (see subsec- 
tion [63]). 

CP [%] 




88 145 231 349 



Fig. 8. — Contributions of different effects to the CP fraction as 
a function of frequency for our best bet model with BH spin a» = 
0.5. Shown are observations (blue error bars), the best bet model 
(solid red line), the same dynamical model computed with zero V 
emissivity {ev = 0) in radiative transfer so that CP is produced by 
Faraday conversion (dot-dashed orange), the same model with zero 
Faraday conversion {pq = 0) (short-dashed brown), and the same 
model with zero Faraday rotation {py = 0) (long-dashed green). 
Emissivity in circular V mode contributes little to the observed 
CP, which is mainly due to Faraday conversion. 



There are several radiative transfer effects that con- 
tribute similarly to the polarized fiuxes. Let us consider 
the production of circular polarization in the fiow. Fig- 
ure [8] shows the consequence of switching off each phys- 
ical effect for our best bet model with spin a* = 0.5. 
The solid red curve is the result with all physics on. The 
dot-dashed orange line below is for zero circular emis- 
sivity having ey = 0. The brown dashed line corre- 



sponds to zero Faraday conversion (pg = 0). Switching 
off ev emissivity leads to a minor correction, whereas 
setting Faraday conversion to zero results in CP of the 
opposite sign with several times smaller absolute value. 
Most of the CP in this model is produced by Faraday 
conversion. It would be incorrect, however, to think 
that the simple linear to circular conversion explains the 
observed CP. The dashed green line in Figure [8] shows 
the CP fraction, when Faraday rotation is switched off 
{pv — 0). The effect of Faraday rotation is insignificant 
at 1/ > 350 GHz, but the rotation of the plane of linear 
polarization simultaneous with conversion between lin- 
ear and circular polarizations produces a unique effect at 
lower I/. This is the so-ca lled "rotation- induced conver- 
sion" (jHoman et al.l [20091 ). Sign oscillations of V with 
frequency do not happen when the Faraday rotation is 
on, but they do happen when py = 0. For the best fit 
model it is the rotation-induced Faraday rotation, which 
is responsible for the most of circularly polarized light. 

In Figure [9] we illustrate the influence of Faraday rota- 
tion on LP fraction (left panel) and EVPA (right panel) . 
The solid curves are produced with all physics on for our 
best bet model with spin a* = 0.5. The green dashed 
lines are computed when switching off Faraday rotation 
{pv = 0) . The Faraday rotation is small at high frequen- 
cies and LP curves look similar at i^ > 200 GHz. As the 
rotation of polarization plane is much stronger at low v, 
a significant phase shift accumulates between different 
rays at the low end of the spectrum and cancellations 
of LP become strong at i/ < 150 GHz. This illustrate s 
the effect of Faraday depolarization ([Bower et al.lll999al ). 
In the absence of Faraday rotation, the dependence of 
EVPA on frequency is not constant: the variations of in- 
trinsic emitted EVPA are significant. Thus, the change 
of EVPA with v should not always be ascribed to the ef- 
fect of Faraday rotation. The positive observed slope of 
EVPA with 1/ at high v, acquired due to negative Fara- 
day rotation measure {RM < 0), is comparable to the 
slope of intrinsic emitted EVPA. 

There is an alternate way to test dynamical models 
against observations. The i ntrinsic im age size was re- 
cently measured (Doeleman et al.l 120081 ) with the VLBI 
technique. The measured correlated flux at 230 GHz was 
Fcorr ~ 0.35 Jy at 3.5 GA SMT-JCMT baseline. Similar 
values of correlated flux were observed later by the same 
group ([Fish et al.l 120 11[ ) . We plot this correlated flux 
with 3ct error bar in Figure [T0| and compare it to simu- 
lated correlated fluxes . To simulate the correlated flux we 
foUow lFish et al.[ ([2009[ ) and employ a Gaussian interstel- 
lar scattering ellipse with half-widths at half-maximum 
7.0 X 3.8GA with position angle 170° East of North. The 
correlated fluxes for the best fit models with spin a^ — 0.5 
(darker red lines) and a* = 0.98 (lighter orange lines) 
are shown. We vary the position angle (PA) of the BH 
spin axis, and plot correlated flux curves with the largest 
(upper solid lines) and the smallest (lower dashed lines) 
correlated flux at 3.5GA. Since we do not fit EVPA di- 
rectly, models with different PA have the same x^/dof. 
The size in our best bet model with spin a^ — 0.5 is 
consistent with observations, whereas the best fit model 
with spin a* = 0.98 has larger correlated flux, so that the 
size of the shadow is slightly under-predicted. The sim- 
ulated source size is in crude agrement to the observed 
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Fig. 9. — Contributions of different effeets to the LP fraction (on the left) and EVPA (on the right) as functions of frequency for the 
best bet model with spin a* = 0.5. Shown are observations (blue error bars and thick blue line), the best bet model (solid red line), and 
the same dynamical model computed with zero Faraday rotation (pv' = 0) in radiative transfer (long-dashed green). Beam depolarization 
is weak: if Faraday rotation is absent, then LP stays high at low frequencies. Even when the Faraday rotation is set to zero, the EVPA 
depends on frequency due to varying intrinsic emission EVPA. Faraday rotation in the best bet model is too weak to reproduce EVPA 
observations, so stronger magnetic fields or more magnetic flux near the black hole than in the simulations may be required. 

over single frozen snapshots, e.g. for Ai = 0. When the 
fast light approximation is used instead of the correct si- 
multaneous evolution of photon field and MHD, the mod- 
els with spins a» — 0; 0.9; 0.98 produce almost identical 
best fits with variations A6' < 0.6°, IST^jT^ < 1.5%, and 
AM/M < 5%. However, the models with a* = 0.5; 0.7 
settle to different x^/dof minima with larger changes in 
quantities: AO = 5°, ATe/Te = 10%, AM/M = 10%. 
These variations are still smaller than variations be- 
tween models with different spins. Switching to the 
fast light approximation results in significant changes 
Ax^/dof ~ 1 between the best fit models for the same 
spins, which emphasizes the importance of precise radi- 
ation transfer calculations. 

6.1. Model Parameters 

We now discuss the estimated parameters obtained 
for the best fit models. The best bet model with spin 
a* = 0.5 has inclination angle 6 = 74.5°, mean accre- 
tion rate M = 4.6 x lO^*Af0year^^, ratio of temper- 
atures Tp/T^ = 20.1 at r = 6M, which gives T^ = 
3.1 • 10^° K at r = 6M in the equatorial plane. The 
best fit models with other spins give the inclination an- 
gles: e = 42°64.5°, 53.5°, 57.2° at a* = 0; 0.7; 0.9; 0.98, 
respectively. Thus, the inclination angle for the 5 mod- 
els lies within 6 = 42° — 75°. Our modeling favors neither 
edge-on nor face-on orientations. The electron temper- 
ature Te at r = 6M is surprisingly uniform over a set 
of best fit models. All 5 best fit models with spins from 
a* = to a, = 0.98 presented in Table [5] have electron 
temperature within the tight range 



baseline [GA] 

Fig. 10. — Correlated flux as a function of baseline at 230 GHz 
normalized to the averaged observed flux at 2.82 Jy for the best 
fit models with spin a» = 0.5 (darker red lines) and a» = 0.98 
(lighter orange lines). The upper solid lines show the smallest size 
(largest correlated flux) over all position angles of BH spin axis, 
and the lower dashed lines show the largest size (smallest correlated 
fl ux) over all posit i on a ngles. An observational results presented 
in'Doelem an et al.l II200S) with Scr error bars at baseline 3.5 GA is 
depicted as a vertical black bar for comparison. The size in our 
best bet model with spin a* = 0.5 is consistent with observations, 
whereas the best fit model with spin a, = 0.98 has larger correlated 
flux, so that the size of the shadow is slightly under-predicted. 

Table [5] summarizes the properties of several best fit 
models. Rows 1 — 5 show the model parameters for best 
fits with spins from a* = to a* = 0.98. The simulated 
spectra are computed every 300M from t = 14, OOOM 
till t = 20, OOOM for At = 60M. Then A^poriod = 21 
spectra are averaged to compare to observations. The 
rows 6—11 show the model parameters for models with 
spin a* = 0.5 for spectra averaged over shorter peri- 
ods. That is, A'pcriod = 21 spectra are computed from 
t = 14, OOOM till t == 15, OOOM for the 1-st short pe- 
riod, while the second short period covers the time in- 
terval from t = 15, OOOA'f tih t = 16, OOOAf , etc. When 
comparing the best fit models with spin a, = 0.5 com- 
puted over different simulation periods, we find varia- 
tions in inclination angle A9 = 3° from the mean, the 
electron temperature AT^/T^ = 10%, and the accretion 
rate AM/M = 30%. The spin position angle varies by 
as much as APA = 30°. 

The last 5 rows in Table [5] show the model parame- 
ters for best fits within the "fast light" approximation. 
In this approximation, simulated spectra are computed 



Te = (3.0 - 4.2; 



X 10i"K. 



(30) 



The accretion rate depends strongly on spin. The model 
with spin a* = has an accretion rate M = 7.0 x 
lO^MQyear"^, which is 5 times larger than the accretion 
rate M ~ 1.4 x lO^M0year~^ for the model with spin 
a* = 0.9. Higher spin values give lower accretion rates. 
A natural outcome of fitting pol arized spectrum is th e PA 
of the BH spin axis. Similar to i Huang et al.l ()2009af ). we 
rely on the observed intrinsic EVPA« 111.5° at 230 GHz 
and EVPAw 146.9° at 349 GHz (see §[2]). For the model 
to fit the difference in EVPA, we add a Faraday rotation 
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screen far froni the BH with constant rotation measure 
(RM). Then we compute the required RM and the intrin- 
sic PA to fit the simulated EVPAs at 230 and 349 GHz. 
The best bet model with a* = 0.5 gives PA = 115.3° 
East of North, whereas the next best fit model with spin 
a* = 0.98 requires PA = 120.3°. However, PA is different 
by 90° between the models with spin a^ — and a* — 0.7, 
which indicates that PA can lie within a wide range. In 
sum, some parameters, such as Tg, are estimated to be in 
narrow ranges, while only order of magnitude estimates 
are available for other parameters, such as M . 

With the estimated orientation of the BH spin axis, we 
can plot an image of average radiation intensity from near 
the event horizon. Figure [TT] shows images of total inten- 
sity /i, for the best bet model with spin a* = 0.5 (upper 
left panel), the best fit for spin a* = 0.98 (lower left 
panel), and LP intensity and CP intensity plots for the 
best bet model with a* = 0.5 (upper right and lower right 
panels, correspondingly). The LP average intensity plot 
was made by averaging U and Q intensities separately 
and then finding the total LP fraction and EVPA. Blue 
(predominant) color on the CP plot depicts the regions 
with negative CP intensity and red (scarce) color depicts 
the regions with positive CP intensity. The total V flux 
from this solution is negative {V < 0). The streamlines 
on the LP plot are aligned with EVPA direction at each 
point. The spin axis is rotated by PA = 115.3° East of 
North for the best bet model with spin a, — 0.5 and by 
PA = 120.3° for the best fit model with spin a* = 0.98. 
The spin axis is inclined at 9 to the line of sight, so that 
the either right (West) or left (East) portions of the flow 
are closer to the observer. The color schemes on all plots 
are nonlinear with corresponding calibration bars plot- 
ted on the sides. The numbers at the top of calibration 
bars denote normalizations. 

7. DISCUSSION AND CONCLUSIONS 

Let us compare our results with estimates of Sgr A* 
accretion flow and BH parameters made by other re- 
searchers. 

Two separate searches for spin based on GRMHD 
numerical simulations have been reported so far: 
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and iDexter et all (|2010l ). 
considered the set of spins 
from a* = 0.5 to a* = 0.98 for 2D GRMHD simulations, 
then fitted the X-Ray flux, the 230 GHz flux, and the 
flux slope at 230 GHz. They found at least one model 
for each spin is crudely consistent with the observations 
(see their tab l e 3), a nd their best bet model has a* — 0.9. 
IDexter et all (pOlOh focused on a set of 3D GRMHD, 
then fltted the 230 GHz flux and size estimates, and 
they provided a table of spin probabilities with a* = 0.9 
having the highest P{a). When we fitted the spectrum 
and LP/CP fractions, the model with a* = 0.5 has the 
lowest x^/dof. As for these two groups, we are also un- 
able to provide a statistically significant constraint on a, . 
Other spin est imate s have been based on analytic mod- 
els. iBroderic k et alJ (1200 9. 2010) favor a* = solutions, 
while iHuang et al.l (|2009bl ) favor a* < 0.9 (although they 
do not explore their full model parameter space). 

Another poorly constrained quantity is the mass 
accretion rate. Our estimate M^st = (1-4 — 

7.0) • lO~*M0year~^ is broad. Acceptable mod- 



els in [Moscibrodzka et al.l ()2009[ ) give M from 0.9 • 
10~^M(7)ye ar~^ to 12- 10~^MM year~\ which agrees with 
our range. IDexter et al.l (|2010 D reported 90% confidence 
interval of M for spin a* = 0.9 solutions, while incorpo- 
rating flow size in x^ analysis. Our estimates have some- 
what higher accretion rates than the range M = 5J^2^ x 
lO-^Mgyear-i (90%) in IDexter et all (|2010D . because 
models with lo wer spin naturally ne ed higher M to fit the 
data. Note that IDexter et al.l (|2009l ) found even lower ac- 
cretion rate M{a^ = 0.9) = (1.0 - 2.3) x lO-^Mgyear"! 
when they assumed equality of proton and electron tem- 
peratures {Tp = Tg). 

In addition to spin and accretion rate, we can try to 
estimate the inclination angle and electron tempera- 
ture Te (Te is reported at r = QM in the equatorial 
plane). Our range is ^cst = 4 2° — 75°, which agrees 
with estimate s by o ther groups. IBroderick et al.l (12009.) : 
IDexter et a ll (I2010D repor t ed 9 y 50°. IH uang et al.l 
(|2009al) and IHuang et aD (|2009b[ ) favor shghtly lower 
9 ~ 40° and 45°, respectiv ely, but they have la r ge er - 
ror bars. To estimate Tg, iMoscibrodzka et al.l ( 20091 ) 
and IDexter et al. (2010) use a constant Tp/Te, whereas 
IHuang et al.l (|2009al) and the present work calculated the 
profile of Te- In al l models, T^ i s a sh allow function of 
radius, which made JDexter et all (|2010[ ) estimate a "com- 
mon" Te = (5.4 ± 3.0) X 10^° K (calculated at some dis- 
tance from the BH center). We measure Te at r = 6M, 
and we obtain a narrower range (likely owing to fitting 
of polarized observations) of Te = (3.0 - 4.2) x lO^"^ K. 

One can use two types of observations to estimate 
the BH spin axis position angle: the 230 GHz cor- 
relat ed flux and the EVPA . Usi ng the correl a ted fl ux 
gave IBroderick eraD (|2009| ) and IDexter et all (|2010l ) a 
result of PA = (-70°)-(-20°) = (110°)-(160°). Us- 
ing the EVPA dat a has given slightly different results: 
iMeve r et al] ()2007D predicts the range PA = 60° - 108° 
whereas Huang gets either PA fti 115° (jHuang et al.l 
I2009bl ) or PA w 140° (Hua ng et ani2009ar i. Our values 
of PA are within the ran ge 85° — 171 ° , whi ch is con- 
sistent with predictions in IMever et al.l ()2007f ) and with 
estimates based upon the observed correlated flux. The 
size of the fl ow may depend substa ntially on the lumi- 
nosity state (jBroderick et al.l l2009| ) or the presence of 
non-thermal structures, spiral waves, and other features. 
In some astrophysical sources, PA is directly known from 
spatially resolved jets, and Sgr A* may be one of such 
sources. A t e ntativ e jet feature was revealed in X-rays by 
iMuno et al.l ()2008D in their fig. 8 showing PAjot = 120°. 
This value is close to PA = 115.3° or PA = 120.3° for 
the best fit models with spins a, = 0.5 and a* = 0.98, 
respectively. 

Besides the estimates of accretion rate and flow prop- 
erties based on the inner flow, there exist estimates based 
on the outer flow. iShcherbakov fc Baganofil ()2010[ ) con- 
structed an inflow-outflow model with conduction and 
stellar winds with radiation matching the X-ray surface 
brightness profile observed by Chandra. Their modeo 
had an accretion rate M = 6 ■ lO~*M0year~^ and elec- 
tron temperature Te = 3.6 x 10^° K at r = 6M, which is 



•^ Note that gra vitational radius is defined as rg 
IShcherbakov fc Baganoff (2010' ). 
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images of total and polarized intensities for best-litting models 
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Fig. 11. — Images of polarized intensities for the best fit models: total intensity for spin a* = 0.98 model (lower left); intensities for 
a» = 0.5 model: total intensity (upper left), linear polarized intensity and streamlines along EVPA (upper right), and circular polarized 
intensity (lower right). Distances are in the units of BH mass M. Images are rotated in the picture plane to correspond to the best spin PA: 
PA = 115.3° for the a* = 0.5 model and PA = 120.3° for the a» = 0.98 model. Individual calibration bars are on the sides of corresponding 
plots. The ill-defined polar region does not contribute significantly to the emission. 



consis tent with present results. iShcherbakov fc BaganoffI 
(|201C1| ) constrained the density in the outer-radial flow 
from X-ray observations, while the present work con- 
strains the density in the inner-radial flow from sub-mni 
observations. The density proflle is then found to be 



p (X r 



0.80-0.90. 



(31) 



The density power-law index (3 lie s be tween 
13 = 1.5 for ADAF flow (jNaravan fc Yil [l995l ) and 
13 = 0.5 for the c onvection-dominated accretion 
flow ([Naravan. Igumensh chev. & Abramowicz 2000; 
IQuataert fc Gruzinovll2000( ). However, the modiflcation 
of the power-law index from the steep ADAF proflle 
is likely due to conduction for Sgr A*, not convection 
(jShcherbakov fc BaganofgiMof) . Newer GRMHD sim- 
ulations of radially-extende d disks show a compa rable 
power-law index for density (jMcKinnev et al.ll2012[ ). 

Our dynamical model has limitations and relies on 
several approximations. More c onvergence testing, like 
done in iMcKinnev et alj (|2012l ). is required to ensure 
the 3D GRMHD simulation results are reliable. The 
amount of initial magnetic flux and the fleld geometry 
might have a pronounced effect on simulation results. For 
example, magnetical l y choked accretion flow s (MCAF) 
(jlgumenshchevl 120081 : IMcKinnev et al.l l2012f) may have 
more desirable properties (such as larger Faraday rota- 
tion as discussed related to Figure [9]) for SgrA* com- 
pared to MRI-dominated disks described in the present 



work. The dependence of the estimated accretion flow 
and the BH parameters on the simulation type and the 
initial setup should be carefully explored in future works. 
The polarization is expected to be able to best highlight 
changes in the magnetic field geometry and strength, and 
so our work is an important stepping stone to distinguish 
whether SgrA* is a classical MRI-dominated disk or an 
MCAF. 

The limited dynamical range of our simulations leads 
to another caveat. We flx electron density rig and tem- 
peratures Tp and Tg in the outer flow and extend them 
down to the event horizon. The slopes of these quan- 
tities break at 25 Af radius, where the power law radial 
extrapolation starts. Thus, the density and temperature 
slopes in the inner flow may need to be determined more 
self-consistently. Future simulations will need to cover 
a larger range of radii and p lasma physics effects, such 
as conduction ([Johnson fc Qu ataert 2007t iSharma et al.l 
I2008t IShcherbakov fc Baganoflai2010iV Simulations with 
larger outer radial boundaries that are run for longer will 
also help to fit the Faraday rotation, which happens for 
the present models partially outside of the simulated do- 
main. A proper simulation of the polar region of the 
flow may be important as well. At present, we artifi- 
cially limit the lowest density and highest temperature 
there. If we do not, then numerical artifacts associated 
with excessive numerical dissipation and heating appear 
(similar to those in iMoscibrodzka et al.l 120091 ) . 
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We found the lowest x^/dof for the model with spin 
a* = 0.5, whereas other groups found a* = and 
a* = 0.9 to provide the best fits in their modeling. So, 
there still appears to be no reliable estimate of BH spin 
for SgrA*. One common shortcoming of recent papers is 
the use of thermal electron distribution. If non-thermal 
electrons provide most of the energy for the sub-mm 
peak, then this would invalid ate all prior spin estimates 
(IShcherbakov fc Huang|[20Tll) . 

The radiative transfer we performed has its short- 
comings. The emissivities in our special syn- 
chrotron approximation provide, e.g., 2% agreement 
with exact emissivities (iLeung. Gammie fc Nobld 120111 : 
IShcherbakov fc Huang|l201lD for & = 20 G. 6>r = 1 rad, 
Te = 6.9 • 10^ K, and observed frequency v — 100 GHz. 
Agreement is better for l arger Tg. Non-polarized 
radiative transfer methods (jMoscibrodzka et al.l 120091 : 
iDexter et al.ll20iol) have an intrinsic error that is com- 
parable with our polarized radiative transfer method for 
the same total emissivity £/, but the error is still 1 — 5%. 

There are other unaccounted sources of error. The 
mass of th e BH in the Galactic Center is known to within 
10% (Ghez et al.l[2008: Gillcssen et al. 2009) and the dis- 
tance is known to 5%. We do not expect these uncertain- 
ties to lead to significant changes in our predictions. A 
shift to slightly lower spin may be able to mimic the effect 
of smaller BH or a BH at larger distance. 

An improvement in observations can lead to further in- 
sights on the fiow and black hole parameters. For exam- 
ple, the detailed comparison of fiux, LP, and CP curves 
in Figure [7] shows that the models with different spins 
have discrepancies at frequencies not yet probed by ob- 
servations. In particular, the CP fractions at 88 GHz 
and 690 GHz are different. The EVPA data need im- 
provement as well. EVPA observations are available at 
230 GHz and 349 GHz, but these frequencies are af- 
fected by Faraday rotation. The observations at higher 
frequency, where the Faraday rotation effect is weaker, 
should provide a better estimate of BH spin axis PA. An- 
other important quantity, LP at 88 GHz, has a largely 
unknown value. Its observations are reported in 2 papers. 
Variations in simulated LP(88GHz) are large between 
the best fit models (see Figure [71). Refinement of the 
observed mean LP(88GHz) could potentially help dis- 
criminate between different spins. A measurement of the 
emitting region size or the correlated fiux is also promis- 
ing. Despite the correlated fiux at 230 GHz being mea- 
sured at the SMT-JCMT 3.5GA baseline, the statistics 
of this measurement need to be improved towards being 
comparable with the statistics of total fiux. The corre- 



lated flux observ ations are currently being accumulated 
(|Fish et al.l 1201 It ). The correlated flux at this baseline 
is exponentially sensitive to the physical fiow size. As 
a caveat, the conclusion on image sizes may depend on 
the behavior of matter in the ill-defined polar regions. 
Our models do not exhibit significant emission from high 
latitudes at 230 GHz (see Figure [TT|) or anywhere above 
88 GHz. 

Future work should incorporate rigorous statisti- 
cal analysis, and such analysis should include tem- 
poral information from the observations. The 
time variability properties can be found from the 
simulations and compared to the observed ones . 
In particular, "iet lags" (f^sef-Zadeh et"al] 120081: 



'Maitra, Markoff & Falcke l2009f) a nd quasi -periodic oscil- 
lations (QPOs) (Gcnz el et al.M 2003: Eck art et al.1 120061 : 
[Mivoshi 2010 1) should be i nvestigated using the simula- 
tions (Dol enceet al.ll2012t ). Also, future 3D GRMHD 
simulations will model more radially extended flows, ac- 
count for ADAF/ADIOS type scale-heights of |/i/r| ~ 1, 
capture outflows, and account for the e ffects of accumu- 
lated magnetic flux near the black hole (jMcKinnev et al.l 
|2012| ). Lastly, for the radiative transfer, adding Comp- 
tonization would be one way to test the quiescent X- 

32— o-i within 2-10 keV 



ray luminosity L » 4 • lO'-'^erg 
(IShcherbakov fc BaganofllMol) . 
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APPENDIX 
RADIATIVE TRANSFER CONVERGENCE 

We have devised a novel code for GR polarized radiative transfer. As with any new code, we need to conduct a 
set of convergence tests to ensure it works accurately. First, we need to come up with metrics for assessing accuracy. 
In the present paper we model fluxes at 7 frequencies between 88 GHz and 857 GHz, LP fractions at 3 frequencies 
and CP fractions at 2 frequencies and define x^ as to characterize the goodness of fit. We employ a similar quantity 
xlf /dof to characterize the accuracy of radiative transfer. We define 

XH/doi ^ - g ^^^^, , (Al) 

where Qi^i are simulated polarized fluxes for one set of radiative transfer parameters and Qi^2 are the fluxes for another 
set. The errors (7{Q) are the observed errors of the mean, and the index i runs through all fitted fluxes, LP, and CP 
fractions. When one of the models fits the data exactly, then x|j/dof coincides with x^/dof. We vary the following 
radiative transfer and dynamical model parameters: 

• number of points P/v along North-South axis and along East- West axis in the picture plane, 

• distance from the center P^s measured in horizon radii r^/, where radiative transfer starts, 
dimensionless scale Pfact of the integration region in the picture plane, 
number of simulated spectra iVpcriods for a single model to compute the mean spectrum, 

• time interval Ai of simultaneous propagation of rays and evolution of numerical simulations, 

• extension power-law slope of density profile Prhopoi 

• extension slope of temperature profile Pupo, 

• extension slope of magnetic field profile Pspo- 

Since fiuctuations and differences in x^/dof between different models reach 1, then values xl^/dof < 0.1 are acceptable, 
but, in general, we strive for x|//dof < 0.02. We set constant Pfact, ^ss, -Psnxy for all radiative transfer computations, 
but we cannot check the code accuracy for all models. We check the convergence a posteriori for the best fit model at 
each spin value. 

We find values of parameters by trial-and-error. The resulting set has Pfact = 1, Pss = l-Olr^, Psnxy — HI, 
-^periods = 21, Ai = 120M. The values of Prhopo and Pupo are fixed by extensions to large radii of temperature and 
density in the inner fiow. 

The tests and the values of xl//dof are summarized in Table[3] The second column describes the test. In particular, 
Pfact : 1 -^ 0.8 means that we test convergence of the integration region relative size. We change one parameter at a 
time. Since the power-law slopes Prhopo and Pupo can vary from model to model, we change them in such a way that 
Prhopo is increased by 0.2 and Pupo is decreased by 0.1. We also estimate the infiuence of magnetic field extension 
power-law slope Pspo by making it shallower from (r/M)^^-'^ to (r/A/)^"-^. We chose to test relatively small variations 
APrhonn = 0.2 and APnnn = 0.1, because dens itv and temperature at Tout = 3 • lO^M are known to within a factor of 
several (jBaganoff et al.ll2003l : iShcherbakov fc Baganofi 2010 ). while these variations correspond to changes by factors 
of 7 and 2.5 in density and temperature, respectively, at rout- 

The results of the tests are as follows. The first 11 tests represent variations of radiative transfer parameters and 
last 3 tests explore the variations of power-law extension slopes. Tests 1 — 4 produce small xff /dof , so that Pn can 
be lowered and Pss can be increased. The changes in the integration region scale (Pfact) result in high x|f/dof ~ 0.07, 
as indicated by tests 5 and 6. Low Pfact leads to systematic underproduction of total flux, whereas high Pfact mainly 
leads to different LP fractions. Test 7 results in high xl^/dof « 0.4, so that a small number of simulated spectra (e.g. 
^periods = H) cauuot be justified. Lower values xlf/dof « 0.13 attained in test 8 indicate that A'^periods = 21 periods 
might be acceptable. With tests 9 — 11, we tested variations in the time interval Aj of simultaneous propagation 
of rays and evolution of numerical simulations. It is expected that longer intervals lead to convergence. However, 
switching from Ai = 120A/ to Ai = 180M and switching Ai = 80Af to Ai = 120A'/ both lead to xl^/dof < 0.1. Since 
these values of xl^/dof are acceptable, we implement Ai = 120A/ for radiative transfer runs. As elucidated by test 
11, freezing simulations in time leads to xl//dof w 0.5, which is too high. Thus, conducting radiative transfer over 
frozen simulation snapshots is not acceptable. Changes in extension slopes of density and temperature (tests 12 and 
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TABLE 3 

Values of x|f/dof for radiative transfer convergence tests and sensitivity to model parameters tests for best fit models. 



Number 


Description 


spin a, = 


spin a, = 0.5 


spin a» = 0.7 


spin a» = 0.9 


spin a» = 0.98 


1 


Pn 


75 -^ 111 


0.00081 


0.00138 


0.00101 


0.00047 


0.01175 


2 


Pn 


111^ 161 


0.00017 


0.00072 


0.00018 


0.00007 


0.00084 


3 


Pss 


l.OOSrH -> l.OlrH 


0.00036 


0.00059 


0.00110 


0.00095 


0.00051 


4 


Pss 


l.Olrn ~* l.OSrn 


0.00778 


0.00982 


0.01616 


0.01468 


0.00920 


5 


Pfact : 0.8 ^ 1.0 


0.01358 


0.07278 


0.06905 


0.02893 


0.02686 


6 


Pfact : 1.0 ^ 1.2 


0.00087 


0.05532 


0.07173 


0.02681 


0.03534 


7 


JVpcHods : 11 -> 21 


0.15665 


0.40611 


0.12233 


0.21397 


0.13588 


8 


Afperiods : 21 ^ 41 


0.02474 


0.04505 


0.13244 


0.02684 


0.04834 


9 


interval At; 120M -^ 180M 


0.06987 


0.06851 


0.13549 


0.02948 


0.10979 


10 


interval At; 80M -5> 120M 


0.09095 


0.04103 


0.03094 


0.03881 


0.06246 


11 


interval At; OM -> 120M 


0.09051 


0.35296 


0.53045 


0.02731 


0.07881 


12 


^rhopo ; Q ^ 0+ = 0.2 


0.04493 


0.04057 


0.03134 


0.01587 


0.05241 


13 


Ajpo : Q^Q- = 0.1 


0.01200 


0.02726 


0.00977 


0.01088 


0.04174 


14 


i^po ; -1.0^ -0.8 


0.02401 


1.05214 


0.15156 


0.05486 


0.04941 



13) result in small Xn/'^oi < 0.05. Variations of magnetic field slope (test 14) lead to large xlf/dof « 1, which means 
the modifications of b extensions will change the best fits. Extensions as shallow as |b| ex (r/M)^°-^ may provide 
better fits to Faraday rotation measure and should be carefully explored. Various extensions of the fluid velocity lead 
to practically the same polarized intensities and are not included in tests. 



